import numpy as np
import matplotlib.pyplot as plt

np.random.seed(50) # Pour assurer la reproductibilité des résultats

Dans ce notebook, on va s'intéresser à l'estimation de la moyenne d'une variable aléatoire donnée. En effet, il s'agit d'un cas particulier d'algorithme d'approximation stochastique.

Soit $X_1,...,X_n$ une suite de variable aléatoire i.i.d de même loi. Si de plus, $\mathbb{E}(X_1)=\mu<\infty,$ la loi des grands nombres nous dit que pour $N$ assez grand, $$\overline{X}_n = \frac{1}{n} \sum_{k=1}^{n} X_k \xrightarrow[\ n \to +\infty \ ]{\ \text{p.s.}\ } \mathbb{E}(X_1).$$

Considérons la suite $\overline{X}_{n+1}$. On a $\forall n\in \mathbb{N},$

$ \begin{align} \overline{X}_{n+1} &= \frac{1}{n+1} \sum_{k=1}^{n+1} X_{k} \\ &= \frac{n}{n+1}\overline{X}_{n}+ \frac{1}{n+1}X_{n+1} \\ &=\overline{X}_n-\frac{1}{n+1}(\overline{X}_n-X_{n+1}). \end{align}$

Dans le cadre d'un algorithme d'approximation stochastique, cela correspond à la suite définie par: $$ \begin{equation} \forall n\geq0, y_{n+1}=y_n-\gamma_{n+1}(y_n-X_{n+1}), \tag{1} \end{equation}$$ où $y_{n+1}=\overline{X}_{n+1}$, $y_n=\overline{X}_{n}$, $\gamma_{n+1}=\frac{1}{n+1}$ et $\mathbb{E}(X_1)=\mu<\infty$.

Si on considère la fonction $$ \begin{align} H:\mathbb{R}^2 &\longrightarrow\mathbb{R} \\ (y,X) &\longmapsto (y-X)^2, \end{align} $$ avec $X$ une suite de variable aléatoire i.i.d de même loi que $(X_n)_{n\geq0}$.

On a alors, $$ \begin{align} h(y) &=\mathbb{E}[H(y,X)]\\ &=\mathbb{E}[y^2-2yX+X^2]\\ &=y^2-2y\mathbb{E}[X]+\mathbb{E}[X^2]. \end{align} $$

La fonction $h$ est donc une foction définie de $\mathbb{R} \mapsto \mathbb{R}$. Elle est continue et dérivable sur $\mathbb{R}$, car il s'agit d'une fonction polynomiale sur $\mathbb{R}$.

En étudiant la fonction $h$, on déduit qu'elle admet un unique minimum et celui-ci est $y^*=\mathbb{E}[X]$.

Approximation Stochastique

Nous voulons maintenant retrouver une approximation de $y^*$ par l'algorithme $(1)$.

En initialisant $y$ à $0$, l'algorithme $(1)$ devient: $$ \begin{equation} \forall n\geq1, y_{n}=\sum_{k=1}^{n}(\prod_{i=k+1}^{n}(1-\gamma_i))\gamma_kX_{k}. \tag{2} \end{equation}$$

D'après le cours, $y_n$ converge vers $y^*$ si et seulement si $\sum_{k\geq1}\gamma_k\xrightarrow[\ n \to +\infty \ ]{\ \text{p.s.}\ }+\infty$ et $\sum_{k\geq1}(\gamma_k)^2<\infty, \text{p.s.}$

Implémentation

On s'intéresse maintenant au cas où $X\sim\mathcal{N}(\mu,\sigma^2)$.

# 1) Paramètres de simulation
N = 10000  # Nombre total de pas de simulation
n_0 = 50  # Numéro de l'étape pour mettre en avant le caractère aléatoire de Y_n0
mu = 5  # Moyenne de la variable aléatoire Z
sigma = 2  # Écart-type de Z
gamma = lambda n: 1 / (n + 1)  # Fonction pour la taille des pas

# Initialisation des tableaux
Y = np.zeros(N + 1)  # Tableau pour stocker les valeurs de Y_n
Y[0] = 0  # Valeur initiale de Y
X_simu = np.random.normal(mu, sigma, N)  # Échantillons i.i.d. de Z, tirés d'une loi normale

# Boucle d'approximation stochastique
for n in range(N):
    # Mise à jour de Y_n avec la règle d'approximation stochastique
    Y[n + 1] = Y[n] - gamma(n) * (Y[n] - X_simu[n])

# 2) Affichage de la valeur de Y_n0 pour souligner l'effet aléatoire
print(f"Valeur de Y_{n_0} (pour montrer le caractère aléatoire): {Y[n_0]}")

# 3) Tracé de la convergence de Y_n vers la moyenne mu
plt.figure(figsize=(14, 6))

# Graphique de convergence
plt.subplot(1, 3, 1)
plt.plot(range(N + 1), Y, label="Y_n", color="blue")  # Tracé de Y_n en fonction de n
plt.axhline(y=mu, color="red", linestyle="--", label="Moyenne réelle (μ)")  # Ligne horizontale pour la moyenne mu
plt.axvline(x=n_0, color="green", linestyle="--", label=f"n = {n_0}")  # Ligne verticale pour n_0
plt.title("Convergence de $Y_n$ vers la moyenne $\\mu$")
plt.xlabel("n")  # Étiquette de l'axe x
plt.ylabel("$Y_n$")  # Étiquette de l'axe y
plt.legend()  # Légende pour les courbes

# 4) Taux de convergence |Y_n - mu|
convergence_rate = np.abs(Y - mu)  # Calcul du taux de convergence
plt.subplot(1, 3, 2)
plt.plot(range(N + 1), convergence_rate, color="purple")  # Tracé du taux de convergence en fonction de n
plt.title("Taux de convergence : $|Y_n - \\mu|$")
plt.xlabel("n")  # Étiquette de l'axe x
plt.ylabel("$|Y_n - \\mu|$")  # Étiquette de l'axe y
plt.yscale("log")  # Échelle logarithmique pour l'axe y

# 5) Tracé de la séquence de taille de pas γ(n)
pas_n = np.array([gamma(n) for n in range(N + 1)])  # Calcul de la séquence de pas γ(n)
plt.subplot(1, 3, 3)
plt.plot(range(N + 1), pas_n, color="orange")  # Tracé de γ(n) en fonction de n
plt.title("Séquence de taille de pas $\\gamma_n = \\frac{1}{n+1}$")
plt.xlabel("n")  # Étiquette de l'axe x
plt.ylabel("$\\gamma_n$")  # Étiquette de l'axe y
plt.yscale("log")  # Échelle logarithmique pour l'axe y

plt.tight_layout()  # Ajustement de la mise en page pour éviter les chevauchements
plt.show()  # Affichage des graphiques

Valeur de Y_50 (pour montrer le caractère aléatoire): 5.186312635066971

# Nombre d'expériences pour illustrer la variabilité à l'étape n_0
n_exp = 30  # Nombre d'expériences à réaliser
Yn_stok = []  # Liste pour stocker les valeurs de Y_{n_0} de chaque expérience

# Répétition de l'approximation stochastique jusqu'à l'étape n_0 pour chaque expérience
for _ in range(n_exp):
    # Réinitialisation de Y et génération d'un nouvel ensemble d'échantillons Z pour chaque expérience
    Y = np.zeros(N + 1)  # Remise à zéro de Y
    X_simu = np.random.normal(mu, sigma, N)  # Nouveaux échantillons i.i.d. pour chaque expérience

    # Exécution de l'approximation uniquement jusqu'à l'étape n_0
    for n in range(n_0):
        Y[n + 1] = Y[n] - gamma(n) * (Y[n] - X_simu[n])  # Mise à jour de Y selon la règle d'approximation stochastique
    
    # Enregistrement du résultat à l'étape n_0
    Yn_stok.append(Y[n_0])

# Tracé de la distribution de Y_{n_0} après n_0 étapes pour les différentes expériences
plt.figure(figsize=(10, 6))
plt.hist(Yn_stok, bins=10, color="skyblue", edgecolor="black", density=True)  # Histogramme des valeurs de Y_{n_0}
plt.axvline(mu, color="red", linestyle="--", label="Moyenne réelle (μ)")  # Ligne indiquant la moyenne réelle
plt.title(f"Distribution de $Y_{{{n_0}}}$ après {n_0} étapes sur {n_exp} expériences")
plt.xlabel(f"$Y_{{{n_0}}}$")  # Étiquette de l'axe x
plt.ylabel("Fréquence")  # Étiquette de l'axe y
plt.legend()  # Légende
plt.show()  # Affichage de l'histogramme

Choix alternatif: $\gamma_n=\frac1{n}$ et $\gamma_n=\frac1{\sqrt{1+n}}$

# Génération d'échantillons i.i.d. de Xi avec une moyenne connue (mu) et une variance (sigma^2) 
mu = 5.0           # Moyenne réelle de la distribution
sigma = 2.0        # Écart-type de la distribution
n_samples = 10000  # Nombre d'échantillons

# Générer des échantillons i.i.d. d'une distribution normale avec moyenne mu et écart-type sigma
Xi = np.random.normal(loc=mu, scale=sigma, size=n_samples)

# Fonction d'estimation récursive de la moyenne avec une séquence de gain générale
def recursive_mean_estimator(Xi, gamma_func):
    Y = np.zeros(len(Xi))  # Tableau pour stocker les estimations récursives
    Y[0] = Xi[0]           # Initialisation avec le premier échantillon comme point de départ
    
    # Parcourir les échantillons et mettre à jour Y avec la formule récursive
    for k in range(1, len(Xi)):
        gamma_k = gamma_func(k)  # Récupérer le gain pour cette itération
        Y[k] = Y[k-1] + gamma_k * (Xi[k] - Y[k-1])  # Mise à jour récursive
    
    return Y

# Définition de différentes fonctions de gain pour comparaison
def gamma_standard(k):
    return 1 / k   # Choix: 1/k

def gamma_sqrt(k):
    return 1 / np.sqrt(k + 1)  # Choix : 1/sqrt(k+1)

# Calcul des estimations récursives avec les deux séquences de gain
Y_standard = recursive_mean_estimator(Xi, gamma_standard)
Y_sqrt = recursive_mean_estimator(Xi, gamma_sqrt)

# Tracé des résultats pour visualiser la convergence de Y_k vers la vraie moyenne
plt.figure(figsize=(12, 6))
plt.plot(Y_standard, label=r'Gain Standard $\gamma_k = \frac{1}{k}$', color='blue')  # Courbe avec le gain standard
plt.plot(Y_sqrt, label=r'Gain Racine Carrée $\gamma_k = \frac{1}{\sqrt{k+1}}$', color='green')  # Courbe avec le gain racine carrée
plt.axhline(y=mu, color='red', linestyle='--', label=r'Moyenne Réelle $\mu = 5$')  # Ligne pour la moyenne réelle
plt.xlabel('Nombre d\'échantillons (k)')  # Étiquette pour l'axe x
plt.ylabel(r'Estimation de $\mu$')  # Étiquette pour l'axe y
plt.legend()  # Légende pour identifier les courbes
plt.title(r'Estimation Récursive de la Moyenne avec Différents Gains $\gamma_k$')  # Titre du graphique
plt.grid(True)  # Activer la grille pour faciliter la lecture
plt.show()  # Afficher le graphique

La figure ci-dessus nous permet de constater que le choix $\gamma_n=\frac1{\sqrt{1+n}}$ n'est pas un choix optimal.